Load the necessary libraries
library(rstanarm) #for fitting models in STAN
library(brms) #for fitting models in STAN
library(coda) #for diagnostics
library(bayesplot) #for diagnostics
library(ggmcmc) #for MCMC diagnostics
library(DHARMa) #for residual diagnostics
library(rstan) #for interfacing with STAN
library(emmeans) #for marginal means etc
library(broom) #for tidying outputs
library(tidybayes) #for more tidying outputs
library(ggeffects) #for partial plots
library(tidyverse) #for data wrangling etc
library(broom.mixed)#for summarising models
library(ggeffects) #for partial effects plots
theme_set(theme_grey()) #put the default ggplot theme back
Here is a modified example from Peake and Quinn (1993). Peake and Quinn (1993) investigated the relationship between the number of individuals of invertebrates living in amongst clumps of mussels on a rocky intertidal shore and the area of those mussel clumps.
Format of peakquinn.csv data files
| AREA | INDIV |
|---|---|
| 516.00 | 18 |
| 469.06 | 60 |
| 462.25 | 57 |
| 938.60 | 100 |
| 1357.15 | 48 |
| ... | ... |
| AREA | Area of mussel clump mm2 - Predictor variable |
| INDIV | Number of individuals found within clump - Response variable |
The aim of the analysis is to investigate the relationship between mussel clump area and the number of non-mussel invertebrate individuals supported in the mussel clump.
peake = read_csv('../public/data/peakquinn.csv', trim_ws=TRUE)
glimpse(peake)
## Rows: 25
## Columns: 2
## $ AREA <dbl> 516.00, 469.06, 462.25, 938.60, 1357.15, 1773.66, 1686.01, 1786…
## $ INDIV <dbl> 18, 60, 57, 100, 48, 118, 148, 214, 225, 283, 380, 278, 338, 27…
When exploring these data as part of a frequentist analysis, exploratory data analysis revealed that the both the response (counds of individuals) and predictor (mussel clump area) were skewed and the relationship between raw counds and mussel clump area was not linear. Furthermore, there was strong evidence of a relationship between mean and variance. Normalising both reponse and predictor addressed these issues. However, rather than log transform the response, it was considered more appropriate to model against a distribution that used a logarithmic link function.
The individual observations here (\(y_i\)) are the observed number of (non mussel individuals found in mussel clump \(i\). As a count, these might be expected to follow a Poisson (or perhaps negative binomial) distribution. In the case of a negative binomial, the observed count for any given mussel clump area are expected to be drawn from a negative binomial distribution with a mean of \(\lambda_i\). All the negative binomial distributions are expected to share the same degree of dispersion (\(\theta\)) - that is, the degree to which the inhabitants of mussell clumps aggregate together (or any other reason for overdispersion) is independent of mussel clump area and can be estimated as a constant across all populations.
The natural log of the expected values (\(\lambda_i\)) is modelled against a linear predictor that includes an intercept (\(\beta_0\)) and a slope (\(\beta_1\)) associated with the natural log of mussel area.
The priors on \(\beta_0\) and \(\beta_1\) should be on the natura log scale (since this will be the scale of the parameters). As starting points, we will consider the following priors:
Model formula: \[ \begin{align} y_i &\sim{} \mathcal{NB}(\lambda_i, \theta)\\ ln(\lambda_i) &= \beta_0 + \beta_1 ln(x_i)\\ \beta_0 & \sim\mathcal{N}(0,5)\\ \beta_1 & \sim\mathcal{N}(0,2)\\ \theta &\sim{} \mathcal{Exp}(1) \end{align} \]
summary(glm(INDIV ~ log(AREA), data=peake, family=poisson()))
##
## Call:
## glm(formula = INDIV ~ log(AREA), family = poisson(), data = peake)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -16.717 -3.614 -1.679 1.808 21.744
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.013528 0.091147 0.148 0.882
## log(AREA) 0.691628 0.009866 70.100 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 7494.5 on 24 degrees of freedom
## Residual deviance: 1547.8 on 23 degrees of freedom
## AIC: 1738.9
##
## Number of Fisher Scoring iterations: 4
summary(MASS::glm.nb(INDIV ~ log(AREA), data=peake))
##
## Call:
## MASS::glm.nb(formula = INDIV ~ log(AREA), data = peake, init.theta = 7.369263003,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.29133 -0.59599 -0.06927 0.48929 1.64734
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.16452 0.53387 -2.181 0.0292 *
## log(AREA) 0.82469 0.06295 13.102 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(7.3693) family taken to be 1)
##
## Null deviance: 161.076 on 24 degrees of freedom
## Residual deviance: 25.941 on 23 degrees of freedom
## AIC: 312.16
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 7.37
## Std. Err.: 2.16
##
## 2 x log-likelihood: -306.16
In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.
peake.rstanarm = stan_glm(INDIV ~ log(AREA), data=peake,
family=poisson(),
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0)
prior_summary(peake.rstanarm)
## Priors for model 'peake.rstanarm'
## ------
## Intercept (after predictors centered)
## ~ normal(location = 0, scale = 2.5)
##
## Coefficients
## Specified prior:
## ~ normal(location = 0, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 0, scale = 2)
## ------
## See help('prior_summary.stanreg') for more details
This tells us:
for the intercept, it is using a normal prior with a mean of 0 and a standard deviation of 2.5. These are the defaults for all non-Gaussian intercepts.
for the coefficients (in this case, just the slope), the default prior is a normal prior centered around 0 with a standard deviation of 2.5. For Poisson models, this is then adjusted for the scale of the data by dividing the 2.5 by the standard deviation of the (in this case log) predictor (then rounded).
2.5/sd(log(peake$AREA))
## [1] 2.025039
One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging preditions would ensure that the priors are unlikely to influence the actual preditions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of paramter space that are not supported by the actual data, the sampler can become unstable and have difficulty.
We can draw from the prior predictive distribution instead of conditioning on the response, by updating the model and indicating prior_PD=TRUE. After refittin the model in this way, we can plot the predictions to gain insights into the range of predictions supported by the priors alone.
peake.rstanarm1 <- update(peake.rstanarm, prior_PD=TRUE)
ggemmeans(peake.rstanarm1, ~AREA) %>% plot(add.data=TRUE)
Conclusions:
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):
I will also overlay the raw data for comparison.
peake.rstanarm2= stan_glm(INDIV ~ log(AREA), data=peake,
family=poisson(),
prior_intercept = normal(0, 5, autoscale=FALSE),
prior = normal(0, 2, autoscale=FALSE),
prior_PD=TRUE,
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0
)
ggemmeans(peake.rstanarm2, ~AREA) %>%
plot(add.data=TRUE)
Now lets refit, conditioning on the data.
peake.rstanarm3= update(peake.rstanarm2, prior_PD=FALSE)
posterior_vs_prior(peake.rstanarm3, color_by='vs', group_by=TRUE,
facet_args=list(scales='free_y'))
Conclusions:
ggemmeans(peake.rstanarm3, ~AREA) %>% plot(add.data=TRUE)
In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.
Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.
peake.brm = brm(bf(INDIV ~ log(AREA), family=poisson()),
data=peake,
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0)
prior_summary(peake.brm)
## prior class coef group resp dpar nlpar bound source
## (flat) b default
## (flat) b logAREA (vectorized)
## student_t(3, 5.8, 2.5) Intercept default
This tells us:
## [1] 5.823046
mad(log(peake$INDIV))
## [1] 1.025465
for the beta coefficients (in this case, just the slope), the default prior is a inproper flat prior. A flat prior essentially means that any value between negative infinity and positive infinity are equally likely. Whilst this might seem reckless, in practice, it seems to work reasonably well for non-intercept beta parameters.
there is no sigma parameter in a binomial model and therefore there are no additional priors.
One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging preditions would ensure that the priors are unlikely to influence the actual preditions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of paramter space that are not supported by the actual data, the sampler can become unstable and have difficulty.
In brms, we can inform the sampler to draw from the prior predictive distribution instead of conditioning on the response, by running the model with the sample_prior='only' argument. Unfortunately, this cannot be applied when there are flat priors (since the posteriors will necessarily extend to negative and positive infinity). Therefore, in order to use this useful routine, we need to make sure that we have defined a proper prior for all parameters.
peake.brm1 = brm(bf(INDIV ~ log(AREA), family=poisson()),
data=peake,
prior=c(
prior(normal(0, 2), class='b')),
sample_prior = 'only',
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0)
ggemmeans(peake.brm1, ~AREA) %>% plot(add.data=TRUE)
conditional_effects(peake.brm1) %>% plot(points=TRUE)
Conclusions:
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):
I will also overlay the raw data for comparison.
peake.brm2 = brm(bf(INDIV ~ log(AREA), family=poisson()),
data=peake,
prior=c(
prior(normal(0, 5), class='Intercept'),
prior(normal(0, 2), class='b')
),
sample_prior = 'only',
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0)
ggemmeans(peake.brm2, ~AREA) %>%
plot(add.data=TRUE)
peake.brm3 <- update(peake.brm2, sample_prior=TRUE, refresh=0)
peake.brm3 %>% get_variables()
## [1] "b_Intercept" "b_logAREA" "prior_Intercept" "prior_b"
## [5] "lp__" "accept_stat__" "stepsize__" "treedepth__"
## [9] "n_leapfrog__" "divergent__" "energy__"
peake.brm3 %>%
posterior_samples %>%
select(-`lp__`) %>%
gather %>%
mutate(Type=ifelse(str_detect(key, 'prior'), 'Prior', 'b'),
Class=ifelse(str_detect(key, 'Intercept'), 'Intercept',
ifelse(str_detect(key, 'b'), 'b', NA))) %>%
ggplot(aes(x=Type, y=value)) +
stat_pointinterval()+
facet_wrap(~Class, scales='free')
standata(peake.brm3)
## $N
## [1] 25
##
## $Y
## [1] 18 60 57 100 48 118 148 214 225 283 380 278 338 274 569
## [16] 509 682 600 978 363 1402 675 1159 1062 632
##
## $K
## [1] 2
##
## $X
## Intercept logAREA
## 1 1 6.246107
## 2 1 6.150731
## 3 1 6.136106
## 4 1 6.844389
## 5 1 7.213142
## 6 1 7.480800
## 7 1 7.430120
## 8 1 7.487896
## 9 1 8.035949
## 10 1 8.289067
## 11 1 8.394989
## 12 1 8.401037
## 13 1 8.513765
## 14 1 8.400853
## 15 1 8.610818
## 16 1 8.919481
## 17 1 8.873303
## 18 1 9.121503
## 19 1 9.223560
## 20 1 9.136445
## 21 1 9.527275
## 22 1 9.918414
## 23 1 10.115073
## 24 1 10.208912
## 25 1 10.170373
## attr(,"assign")
## [1] 0 1
##
## $prior_only
## [1] 0
##
## attr(,"class")
## [1] "standata"
stancode(peake.brm3)
## // generated with brms 2.14.0
## functions {
## }
## data {
## int<lower=1> N; // total number of observations
## int Y[N]; // response variable
## int<lower=1> K; // number of population-level effects
## matrix[N, K] X; // population-level design matrix
## int prior_only; // should the likelihood be ignored?
## }
## transformed data {
## int Kc = K - 1;
## matrix[N, Kc] Xc; // centered version of X without an intercept
## vector[Kc] means_X; // column means of X before centering
## for (i in 2:K) {
## means_X[i - 1] = mean(X[, i]);
## Xc[, i - 1] = X[, i] - means_X[i - 1];
## }
## }
## parameters {
## vector[Kc] b; // population-level effects
## real Intercept; // temporary intercept for centered predictors
## }
## transformed parameters {
## }
## model {
## // likelihood including all constants
## if (!prior_only) {
## target += poisson_log_glm_lpmf(Y | Xc, Intercept, b);
## }
## // priors including all constants
## target += normal_lpdf(b | 0, 2);
## target += normal_lpdf(Intercept | 0, 5);
## }
## generated quantities {
## // actual population-level intercept
## real b_Intercept = Intercept - dot_product(means_X, b);
## // additionally draw samples from priors
## real prior_b = normal_rng(0,2);
## real prior_Intercept = normal_rng(0,5);
## }
In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.
There are a wide variety of tests that range from the big picture, overal chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.
The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.
available_mcmc()
## bayesplot MCMC module:
## mcmc_acf
## mcmc_acf_bar
## mcmc_areas
## mcmc_areas_data
## mcmc_areas_ridges
## mcmc_areas_ridges_data
## mcmc_combo
## mcmc_dens
## mcmc_dens_chains
## mcmc_dens_chains_data
## mcmc_dens_overlay
## mcmc_hex
## mcmc_hist
## mcmc_hist_by_chain
## mcmc_intervals
## mcmc_intervals_data
## mcmc_neff
## mcmc_neff_data
## mcmc_neff_hist
## mcmc_nuts_acceptance
## mcmc_nuts_divergence
## mcmc_nuts_energy
## mcmc_nuts_stepsize
## mcmc_nuts_treedepth
## mcmc_pairs
## mcmc_parcoord
## mcmc_parcoord_data
## mcmc_rank_hist
## mcmc_rank_overlay
## mcmc_recover_hist
## mcmc_recover_intervals
## mcmc_recover_scatter
## mcmc_rhat
## mcmc_rhat_data
## mcmc_rhat_hist
## mcmc_scatter
## mcmc_trace
## mcmc_trace_data
## mcmc_trace_highlight
## mcmc_violin
Of these, we will focus on:
plot(peake.rstanarm3, plotfun='mcmc_trace')
The chains appear well mixed and very similar
plot(peake.rstanarm3, 'acf_bar')
There is no evidence of autocorrelation in the MCMC samples
plot(peake.rstanarm3, 'rhat_hist')
All Rhat values are below 1.05, suggesting the chains have converged.
neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
plot(peake.rstanarm3, 'neff_hist')
Ratios all very high.
plot(peake.rstanarm3, 'combo')
plot(peake.rstanarm3, 'violin')
The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
stan_trace(peake.rstanarm3)
The chains appear well mixed and very similar
stan_ac(peake.rstanarm3)
There is no evidence of autocorrelation in the MCMC samples
stan_rhat(peake.rstanarm3)
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
stan_ess(peake.rstanarm3)
Ratios all very high.
stan_dens(peake.rstanarm3, separate_chains = TRUE)
The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
peake.ggs <- ggs(peake.rstanarm3)
ggs_traceplot(peake.ggs)
The chains appear well mixed and very similar
ggs_autocorrelation(peake.ggs)
There is no evidence of autocorrelation in the MCMC samples
ggs_Rhat(peake.ggs)
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
ggs_effective(peake.ggs)
Ratios all very high.
ggs_crosscorrelation(peake.ggs)
ggs_grb(peake.ggs)
The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.
available_mcmc()
## bayesplot MCMC module:
## mcmc_acf
## mcmc_acf_bar
## mcmc_areas
## mcmc_areas_data
## mcmc_areas_ridges
## mcmc_areas_ridges_data
## mcmc_combo
## mcmc_dens
## mcmc_dens_chains
## mcmc_dens_chains_data
## mcmc_dens_overlay
## mcmc_hex
## mcmc_hist
## mcmc_hist_by_chain
## mcmc_intervals
## mcmc_intervals_data
## mcmc_neff
## mcmc_neff_data
## mcmc_neff_hist
## mcmc_nuts_acceptance
## mcmc_nuts_divergence
## mcmc_nuts_energy
## mcmc_nuts_stepsize
## mcmc_nuts_treedepth
## mcmc_pairs
## mcmc_parcoord
## mcmc_parcoord_data
## mcmc_rank_hist
## mcmc_rank_overlay
## mcmc_recover_hist
## mcmc_recover_intervals
## mcmc_recover_scatter
## mcmc_rhat
## mcmc_rhat_data
## mcmc_rhat_hist
## mcmc_scatter
## mcmc_trace
## mcmc_trace_data
## mcmc_trace_highlight
## mcmc_violin
Of these, we will focus on:
mcmc_plot(peake.brm3, type='trace')
The chains appear well mixed and very similar
mcmc_plot(peake.brm3, type='acf_bar')
There is no evidence of autocorrelation in the MCMC samples
mcmc_plot(peake.brm3, type='rhat_hist')
All Rhat values are below 1.05, suggesting the chains have converged.
neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
mcmc_plot(peake.brm3, type='neff_hist')
Ratios all very high.
mcmc_plot(peake.brm3, type='combo')
mcmc_plot(peake.brm3, type='violin')
The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
stan_trace(peake.brm3$fit)
The chains appear well mixed and very similar
stan_ac(peake.brm3$fit)
There is no evidence of autocorrelation in the MCMC samples
stan_rhat(peake.brm3$fit)
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
stan_ess(peake.brm3$fit)
Ratios all very high.
stan_dens(peake.brm3$fit, separate_chains = TRUE)
The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
peake.ggs <- ggs(peake.brm3)
ggs_traceplot(peake.ggs)
The chains appear well mixed and very similar
ggs_autocorrelation(peake.ggs)
There is no evidence of autocorrelation in the MCMC samples
ggs_Rhat(peake.ggs)
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
ggs_effective(peake.ggs)
Ratios all very high.
ggs_crosscorrelation(peake.ggs)
ggs_grb(peake.ggs)
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
available_ppc()
## bayesplot PPC module:
## ppc_bars
## ppc_bars_grouped
## ppc_boxplot
## ppc_data
## ppc_dens
## ppc_dens_overlay
## ppc_ecdf_overlay
## ppc_error_binned
## ppc_error_hist
## ppc_error_hist_grouped
## ppc_error_scatter
## ppc_error_scatter_avg
## ppc_error_scatter_avg_vs_x
## ppc_freqpoly
## ppc_freqpoly_grouped
## ppc_hist
## ppc_intervals
## ppc_intervals_data
## ppc_intervals_grouped
## ppc_loo_intervals
## ppc_loo_pit
## ppc_loo_pit_overlay
## ppc_loo_pit_qq
## ppc_loo_ribbon
## ppc_ribbon
## ppc_ribbon_data
## ppc_ribbon_grouped
## ppc_rootogram
## ppc_scatter
## ppc_scatter_avg
## ppc_scatter_avg_grouped
## ppc_stat
## ppc_stat_2d
## ppc_stat_freqpoly_grouped
## ppc_stat_grouped
## ppc_violin_grouped
pp_check(peake.rstanarm3, plotfun='dens_overlay')
The model draws appear deviate from the observed data.
pp_check(peake.rstanarm3, plotfun='error_scatter_avg')
The predictive error seems to be related to the predictor - the model performs poorest at higher mussel clump areas.
pp_check(peake.rstanarm3, x=peake$AREA, plotfun='error_scatter_avg_vs_x')
pp_check(peake.rstanarm3, x=peake$AREA, plotfun='intervals')
The modelled preditions seem to underestimate the uncertainty with increasing mussel clump area.
pp_check(peake.rstanarm3, x=peake$AREA, plotfun='ribbon')
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
#library(shinystan)
#launch_shinystan(peake.rstanarm3)
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
preds <- posterior_predict(peake.rstanarm3, nsamples=250, summary=FALSE)
peake.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = peake$INDIV,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE)
plot(peake.resids)
Conclusions:
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
available_ppc()
## bayesplot PPC module:
## ppc_bars
## ppc_bars_grouped
## ppc_boxplot
## ppc_data
## ppc_dens
## ppc_dens_overlay
## ppc_ecdf_overlay
## ppc_error_binned
## ppc_error_hist
## ppc_error_hist_grouped
## ppc_error_scatter
## ppc_error_scatter_avg
## ppc_error_scatter_avg_vs_x
## ppc_freqpoly
## ppc_freqpoly_grouped
## ppc_hist
## ppc_intervals
## ppc_intervals_data
## ppc_intervals_grouped
## ppc_loo_intervals
## ppc_loo_pit
## ppc_loo_pit_overlay
## ppc_loo_pit_qq
## ppc_loo_ribbon
## ppc_ribbon
## ppc_ribbon_data
## ppc_ribbon_grouped
## ppc_rootogram
## ppc_scatter
## ppc_scatter_avg
## ppc_scatter_avg_grouped
## ppc_stat
## ppc_stat_2d
## ppc_stat_freqpoly_grouped
## ppc_stat_grouped
## ppc_violin_grouped
pp_check(peake.brm3, type='dens_overlay')
The model draws appear deviate from the observed data.
pp_check(peake.brm3, type='error_scatter_avg')
The predictive error seems to be related to the predictor - the model performs poorest at higher mussel clump areas.
pp_check(peake.brm3, x='AREA', type='error_scatter_avg_vs_x')
pp_check(peake.brm3, x='AREA', type='intervals')
The modelled preditions seem to underestimate the uncertainty with increasing mussel clump area.
pp_check(peake.brm3, x='AREA', type='ribbon')
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
#library(shinystan)
#launch_shinystan(peake.brm3)
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
preds <- posterior_predict(peake.brm3, nsamples=250, summary=FALSE)
peake.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = peake$INDIV,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE)
plot(peake.resids)
Conclusions:
peake.rstanarm4= stan_glm(INDIV ~ log(AREA), data=peake,
family=neg_binomial_2(),
prior_intercept = normal(0, 5, autoscale=FALSE),
prior = normal(0, 2, autoscale=FALSE),
prior_aux = rstanarm::exponential(rate=1, autoscale=FALSE),
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0
)
posterior_vs_prior(peake.rstanarm4, color_by='vs', group_by=TRUE,
facet_args=list(scales='free_y'))
Conclusions:
ggemmeans(peake.rstanarm4, ~AREA) %>% plot(add.data=TRUE)
plot(peake.rstanarm4, plotfun='mcmc_trace')
plot(peake.rstanarm4, 'acf_bar')
plot(peake.rstanarm4, 'rhat_hist')
plot(peake.rstanarm4, 'neff_hist')
pp_check(peake.rstanarm4, x=peake$AREA, plotfun='error_scatter_avg_vs_x')
pp_check(peake.rstanarm4, x=peake$AREA, plotfun='intervals')
preds <- posterior_predict(peake.rstanarm4, nsamples=250, summary=FALSE)
peake.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = peake$INDIV,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE)
plot(peake.resids)
Different models can be compared via information criterion. LOO (Leave One Out) IC is similar to AIC except that AIC does not consider priors and assumes that the posterior likelihood is multivariate normal. LOO AIC does not and integrates over all uncertainty.
elpd_loo: is the Bayesian LOO estimate of ELPDSE of elpd_loo: is the Monte Carlo standard error is the estimate for the computational accuracy of MCMC and importance sampling used to compute elpd_loop_loo: is the difference between elpd_loo and the estimated from non-cross validated ELPD and therefore is a measure of the relative extra difficulty of predicting new data compared to the observed data. It can also be a meaure of the effective number of parameters:
p_loo will be less than the number of estimated parameters and the total sample size.p_loo will be greater than the number of estimated parameters and the total sample size.Pareto k: is a relative measure of the influence of each observation as well as the accuracy with which the associated leave one out calculations can be estimated.
k < 0.5: the corresponding components can be accurately estimated0.5 < k < 0.7: the accuracy is lower but still acceptablek>0.7: the accuracy is too low and elpd_loo is unreliable. This can also suggest a misspecified model.More information about the above can be gleaned from [https://mc-stan.org/loo/reference/loo-glossary.html].
Start with the Poisson model
(peake.rstanarm3.loo =loo(peake.rstanarm3))
##
## Computed from 2400 by 25 log-likelihood matrix
##
## Estimate SE
## elpd_loo -932.9 298.9
## p_loo 122.2 49.0
## looic 1865.9 597.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 14 56.0% 869
## (0.5, 0.7] (ok) 2 8.0% 81
## (0.7, 1] (bad) 5 20.0% 20
## (1, Inf) (very bad) 4 16.0% 1
## See help('pareto-k-diagnostic') for details.
Conclusions:
p_loo is substantially higher than the number of estimated parametersPareto k values identified as bad or even very badNow for the Negative Binomial model
(peake.rstanarm4.loo =loo(peake.rstanarm4))
##
## Computed from 2400 by 25 log-likelihood matrix
##
## Estimate SE
## elpd_loo -156.7 5.4
## p_loo 1.9 0.5
## looic 313.4 10.7
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
Conclusions:
p_loo is lower than the number of estimated parametersPareto k valuesWe can also compare the models.
The difference in Expected Log Pointwise Probability Density (elpd_dff) computes the difference in elpd_loo of two models (expressed relative to the higher model). The difference in elpd_diff will be negative if the expected out-of-sample predictive accuracy of the first model is higher. If the difference is be positive then the second model is preferred.
Note, the above assumes that both models are valid. In the current case, we know that the Poisson model was overdispersed and thus it should not be a genuine candidate. Nevertheless, we will compare the Poisson to the Negative Binomial for illustrative purposes.
loo_compare(peake.rstanarm3.loo, peake.rstanarm4.loo)
## elpd_diff se_diff
## peake.rstanarm4 0.0 0.0
## peake.rstanarm3 -776.2 295.4
Conclusions:
elpd_diff is negative, indicating that the first model (Negative Binomial) is better.peake.brm4 = brm(bf(INDIV ~ log(AREA), family=negbinomial()),
data=peake,
prior=c(
prior(normal(0, 5), class='Intercept'),
prior(normal(0, 2), class='b'),
prior(gamma(0.01, 0.01), class='shape')
),
sample_prior=TRUE,
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0)
peake.brm4 %>%
posterior_samples %>%
select(-`lp__`) %>%
gather %>%
mutate(Type=ifelse(str_detect(key, 'prior'), 'Prior', 'b'),
Class=ifelse(str_detect(key, 'Intercept'), 'Intercept',
ifelse(str_detect(key, 'b'), 'b', 'shape'))) %>%
ggplot(aes(x=Type, y=value)) +
stat_pointinterval()+
facet_wrap(~Class, scales='free')
Conclusions:
ggemmeans(peake.brm4, ~AREA) %>% plot(add.data=TRUE)
plot(peake.brm4, type='mcmc_trace')
plot(peake.brm4, type='acf_bar')
plot(peake.brm4, type='rhat_hist')
plot(peake.brm4, type='neff_hist')
pp_check(peake.brm4, x='AREA',type='error_scatter_avg_vs_x')
pp_check(peake.brm4, x='AREA', type='intervals')
preds <- posterior_predict(peake.brm4, nsamples=250, summary=FALSE)
peake.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = peake$INDIV,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE)
plot(peake.resids)
Different models can be compared via information criterion. LOO (Leave One Out) IC is similar to AIC except that AIC does not consider priors and assumes that the posterior likelihood is multivariate normal. LOO AIC does not and integrates over all uncertainty.
elpd_loo: is the Bayesian LOO estimate of ELPDSE of elpd_loo: is the Monte Carlo standard error is the estimate for the computational accuracy of MCMC and importance sampling used to compute elpd_loop_loo: is the difference between elpd_loo and the estimated from non-cross validated ELPD and therefore is a measure of the relative extra difficulty of predicting new data compared to the observed data. It can also be a meaure of the effective number of parameters:
p_loo will be less than the number of estimated parameters and the total sample size.p_loo will be greater than the number of estimated parameters and the total sample size.Pareto k: is a relative measure of the influence of each observation as well as the accuracy with which the associated leave one out calculations can be estimated.
k < 0.5: the corresponding components can be accurately estimated0.5 < k < 0.7: the accuracy is lower but still acceptablek>0.7: the accuracy is too low and elpd_loo is unreliable. This can also suggest a misspecified model.More information about the above can be gleaned from [https://mc-stan.org/loo/reference/loo-glossary.html].
Start with the Poisson model
(peake.brm3.loo =loo(peake.brm3))
##
## Computed from 2400 by 25 log-likelihood matrix
##
## Estimate SE
## elpd_loo -928.4 296.9
## p_loo 120.3 48.2
## looic 1856.7 593.8
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 15 60.0% 522
## (0.5, 0.7] (ok) 3 12.0% 73
## (0.7, 1] (bad) 4 16.0% 24
## (1, Inf) (very bad) 3 12.0% 1
## See help('pareto-k-diagnostic') for details.
Conclusions:
p_loo is substantially higher than the number of estimated parametersPareto k values identified as bad or even very badNow for the Negative Binomial model
(peake.brm4.loo =loo(peake.brm4))
##
## Computed from 2400 by 25 log-likelihood matrix
##
## Estimate SE
## elpd_loo -156.4 5.8
## p_loo 3.0 0.9
## looic 312.7 11.5
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 24 96.0% 720
## (0.5, 0.7] (ok) 1 4.0% 1694
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
Conclusions:
p_loo is lower than the number of estimated parametersPareto k valuesWe can also compare the models.
The difference in Expected Log Pointwise Probability Density (elpd_dff) computes the difference in elpd_loo of two models (expressed relative to the higher model). The difference in elpd_diff will be negative if the expected out-of-sample predictive accuracy of the first model is higher. If the difference is be positive then the second model is preferred.
Note, the above assumes that both models are valid. In the current case, we know that the Poisson model was overdispersed and thus it should not be a genuine candidate. Nevertheless, we will compare the Poisson to the Negative Binomial for illustrative purposes.
loo_compare(peake.brm3.loo, peake.brm4.loo)
## elpd_diff se_diff
## peake.brm4 0.0 0.0
## peake.brm3 -772.0 292.6
Conclusions:
elpd_diff is negative, indicating that the first model (Negative Binomial) is better.peake.rstanarm4 %>% ggpredict() %>% plot(add.data=TRUE)
## $AREA
peake.rstanarm4 %>% ggemmeans(~AREA, type='fixed', transform='response') %>% plot(add.data=TRUE)
peake.rstanarm4 %>% fitted_draws(newdata=peake) %>%
median_hdci() %>%
ggplot(aes(x=AREA, y=.value)) +
geom_ribbon(aes(ymin=.lower, ymax=.upper), fill='blue', alpha=0.3) +
geom_line() +
geom_point(data=peake, aes(y=INDIV, x=AREA))
peake.brm4 %>% conditional_effects()
peake.brm4 %>% conditional_effects(spaghetti=TRUE,nsamples=200)
peake.brm4 %>% ggpredict() %>% plot(add.data=TRUE)
## $AREA
peake.brm4 %>% ggemmeans(~AREA) %>% plot(add.data=TRUE)
peake.brm4 %>% fitted_draws(newdata=peake) %>%
median_hdci() %>%
ggplot(aes(x=AREA, y=.value)) +
geom_ribbon(aes(ymin=.lower, ymax=.upper), fill='blue', alpha=0.3) +
geom_line() +
geom_point(data=peake, aes(y=INDIV, x=AREA))
rstanarm captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.
The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
summary(peake.rstanarm4)
##
## Model Info:
## function: stan_glm
## family: neg_binomial_2 [log]
## formula: INDIV ~ log(AREA)
## algorithm: sampling
## sample: 2400 (posterior sample size)
## priors: see help('prior_summary')
## observations: 25
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) -1.2 0.7 -2.1 -1.2 -0.2
## log(AREA) 0.8 0.1 0.7 0.8 0.9
## reciprocal_dispersion 4.7 1.3 3.1 4.5 6.4
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 479.7 89.0 377.1 471.2 596.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 2458
## log(AREA) 0.0 1.0 2477
## reciprocal_dispersion 0.0 1.0 2343
## mean_PPD 1.7 1.0 2624
## log-posterior 0.0 1.0 2348
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
Conclusions:
tidyMCMC(peake.rstanarm4$stanfit, estimate.method='median', conf.int=TRUE, conf.method='HPDinterval', rhat=TRUE, ess=TRUE)
Conclusions:
Due to the presence of a log transform in the predictor, it is better to use the regex version.
peake.rstanarm4 %>% get_variables()
## [1] "(Intercept)" "log(AREA)" "reciprocal_dispersion"
## [4] "accept_stat__" "stepsize__" "treedepth__"
## [7] "n_leapfrog__" "divergent__" "energy__"
peake.draw <- peake.rstanarm4 %>% gather_draws(`.Intercept.*|.*AREA.*`, regex=TRUE)
peake.draw
We can then summarise this
peake.draw %>% median_hdci
peake.rstanarm4 %>%
gather_draws(`.Intercept.*|.*AREA.*`, regex=TRUE) %>%
ggplot() +
stat_halfeye(aes(x=.value, y=.variable)) +
facet_wrap(~.variable, scales='free')
We could alternatively express the parameters on the response scale.
peake.rstanarm4 %>%
gather_draws(`.Intercept.*|.*AREA.*`, regex=TRUE) %>%
group_by(.variable) %>%
mutate(.value=exp(.value)) %>%
median_hdci
Conclusions:
peake.rstanarm4 %>% plot(plotfun='mcmc_intervals')
This is purely a graphical depiction on the posteriors.
peake.rstanarm4 %>% tidy_draws()
peake.rstanarm4 %>% spread_draws(`.Intercept.*|.*AREA.*`, regex=TRUE)
peake.rstanarm4 %>% posterior_samples() %>% as_tibble()
Unfortunately, \(R^2\) calculations for models other than Gaussian and Binomial have not yet been implemented for rstanarm models yet.
#peake.rstanarm4 %>% bayes_R2() %>% median_hdci
brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.
The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
summary(peake.brm4)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: INDIV ~ log(AREA)
## Data: peake (Number of observations: 25)
## Samples: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup samples = 2400
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.15 0.61 -2.39 0.03 1.00 2165 2033
## logAREA 0.82 0.07 0.68 0.97 1.00 2162 2198
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 6.85 2.12 3.53 11.43 1.00 2217 2315
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Conclusions:
tidyMCMC(peake.brm4$fit, estimate.method='median', conf.int=TRUE, conf.method='HPDinterval', rhat=TRUE, ess=TRUE)
Conclusions:
Due to the presence of a log transform in the predictor, it is better to use the regex version.
peake.brm4 %>% get_variables()
## [1] "b_Intercept" "b_logAREA" "shape" "prior_Intercept"
## [5] "prior_b" "prior_shape" "lp__" "accept_stat__"
## [9] "stepsize__" "treedepth__" "n_leapfrog__" "divergent__"
## [13] "energy__"
peake.draw <- peake.brm4 %>% gather_draws(`b_.*`, regex=TRUE)
peake.draw
We can then summarise this
peake.draw %>% median_hdci
peake.brm4 %>%
gather_draws(`b_.*`, regex=TRUE) %>%
ggplot() +
stat_halfeye(aes(x=.value, y=.variable)) +
facet_wrap(~.variable, scales='free')
We could alternatively express the parameters on the response scale.
peake.brm4 %>%
gather_draws(`.Intercept.*|.*AREA.*`, regex=TRUE) %>%
group_by(.variable) %>%
mutate(.value=exp(.value)) %>%
median_hdci
Conclusions:
peake.brm4 %>% plot(plotfun='mcmc_intervals')
This is purely a graphical depiction on the posteriors.
peake.brm4 %>% tidy_draws()
peake.brm4 %>% spread_draws(`b_.*`, regex=TRUE)
peake.brm4 %>% posterior_samples() %>% as_tibble()
peake.brm4 %>% bayes_R2(summary=FALSE) %>% median_hdci
Since we have the entire posterior, we are able to make probablity statements. We simply count up the number of MCMC sample draws that satisfy a condition (e.g represent a slope greater than 0) and then divide by the total number of MCMC samples.
For this exersize, we will explore the following:
Unfortunately, the inline log-transformation of the predictor interfers with hypothesis()’s ability to find the slope. We can get around this by renaming before calling hypothesis().
peake.rstanarm4 %>% as.data.frame() %>% rename(lAREA=`log(AREA)`) %>% hypothesis('lAREA>0')
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (lAREA) > 0 0.82 0.09 0.69 0.96 Inf 1 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Conclusions:
Alternatively…
peake.rstanarm4 %>% tidy_draws() %>% summarise(P=sum(`log(AREA)`>0)/n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 1
Conclusions:
newdata = list(AREA=c(5000, 10000))
peake.rstanarm4 %>% emmeans(~AREA, at=newdata) %>% pairs()
## contrast estimate lower.HPD upper.HPD
## 5000 - 10000 -0.571 -0.691 -0.459
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95
Conclusions:
If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.
peake.mcmc <- peake.rstanarm4 %>% emmeans(~AREA, at=newdata) %>%
tidy_draws() %>%
rename_with(~str_replace(., 'AREA ', 'p')) %>%
mutate(Eff=p10000 - p5000,
PEff=100*Eff/p5000)
peake.mcmc %>% head
Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).
peake.mcmc %>% tidyMCMC(estimate.method='median',
conf.int=TRUE, conf.method='HPDinterval')
Alternatively, we could use median_hdci
peake.mcmc %>% median_hdci(PEff)
Conclusions:
To get the probability that the effect is greater than a 50% increase.
peake.mcmc %>% summarise(P=sum(PEff>50)/n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 0
Conclusions:
Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).
peake.mcmc %>% hypothesis('PEff>50')
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (PEff)-(50) > 0 -40.26 1 -41.92 -38.66 0 0
## Star
## 1
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Unfortunately, the inline log-transformation of the predictor interfers with hypothesis()’s ability to find the slope.
peake.brm4 %>% hypothesis('logAREA>0')
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (logAREA) > 0 0.82 0.07 0.71 0.94 Inf 1 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Conclusions:
Alternatively…
peake.brm4 %>% tidy_draws() %>% summarise(P=sum(b_logAREA>0)/n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 1
Conclusions:
newdata = list(AREA=c(5000, 10000))
peake.brm4 %>% emmeans(~AREA, at=newdata) %>% pairs()
## contrast estimate lower.HPD upper.HPD
## 5000 - 10000 -0.57 -0.672 -0.476
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95
Conclusions:
If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.
peake.mcmc <- peake.brm4 %>% emmeans(~AREA, at=newdata) %>%
tidy_draws() %>%
rename_with(~str_replace(., 'AREA ', 'p')) %>%
mutate(Eff=p10000 - p5000,
PEff=100*Eff/p5000)
peake.mcmc %>% head
Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).
peake.mcmc %>% tidyMCMC(estimate.method='median',
conf.int=TRUE, conf.method='HPDinterval')
Alternatively, we could use median_hdci
peake.mcmc %>% median_hdci(PEff)
Conclusions:
To get the probability that the effect is greater than a 50% increase.
peake.mcmc %>% summarise(P=sum(PEff>50)/n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 0
Conclusions:
Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).
peake.mcmc %>% hypothesis('PEff>50')
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (PEff)-(50) > 0 -40.26 0.84 -41.63 -38.85 0 0
## Star
## 1
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## Using emmeans
peake.grid = with(peake, list(AREA = seq(min(AREA), max(AREA), len=100)))
newdata = emmeans(peake.rstanarm4, ~AREA, at=peake.grid, type='response') %>% as.data.frame
head(newdata)
ggplot(newdata, aes(y=prob, x=AREA)) +
geom_point(data=peake, aes(y=INDIV)) +
geom_line() +
geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD), fill='blue', alpha=0.3) +
scale_y_continuous('Individuals') +
scale_x_continuous('Mussel clump area') +
theme_classic()
## Using emmeans
peake.grid = with(peake, list(AREA = seq(min(AREA), max(AREA), len=100)))
newdata = emmeans(peake.brm4, ~AREA, at=peake.grid, type='response') %>% as.data.frame
head(newdata)
ggplot(newdata, aes(y=prob, x=AREA)) +
geom_point(data=peake, aes(y=INDIV)) +
geom_line() +
geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD), fill='blue', alpha=0.3) +
scale_y_continuous('Individuals') +
scale_x_continuous('Mussel clump area') +
theme_classic()
Peake, A. J., and G. P. Quinn. 1993. “Temporal Variation in Species-Area Curves for Invertebrates in Clumps of an Intertidal Mussel.” Ecography 16: 269–77.